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Introduction 

The purpose of this research effort is to develop a means to use, and to ultimately 
implement, /ip- version finite elements in the numerical solution of optimal control problems. 
Under NAG-939, the hybrid MACS YM A/FORTRAN code GENCODE was developed 
which utilized /ir version finite elements to successfully approximate solutions to a wide 
class of optimal control problems. In that code the means for improvement of the solution 
was the refinement of the time-discretization mesh. With the extension to /tp- version 
finite elements, the degrees of freedom include both nodal values and extra interior values 
associated with the unknown states, co-states, and controls, the number of which depends 
on the order of the shape functions in each element. For details, see [1]. 

Progress During This Period 


Optimal Control Problems 

A FORTRAN code has been developed using higher-order finite-elements to approx- 
imate solutions to a particular subset of optimal control problems. The cost function to 
be minimized for these problems can contain both an scalar penalty, <j>, on the states, x, 
at the initial time or final time, tf, plus an integral penalty, L, on the states and controls, 
u, in this form: 

J = 4>[x(t 0 ),x(tf),tf] + [ L(x,u) dt (1) 

Jo 

The state rates are governed by differential equations 

x = f(x,u ) x E R n *,u E R nu (2) 

where / is an autonomous function of the state vector x. The boundary conditions can be 
specified at the initial time, the final time, or some combination of both, in the form 

<P[x(0),x(t/),i/] = 0 E R nbc (3) 

This formulation also allows for control inequality constraints of the form: 


g(x , tt) < 0 g E R np n p <n u 


( 4 ) 


which are enforced through use of slack variables, k. 

Adjoining the differential equations, boundary conditions, and control constraints to 
the original cost function by means of Lagrange multipliers A, u, and // respectively yields 
a new cost function J'\ 


J' = <£[x(f 0 ),x(t/),t/] + v T V[x{G),x(tf),tf\ + a T (x - x)|,° + 



+ A T [/(x, it) - x] 4- n T \g(x, u) -F k 2 ] } dt 


( 5 ) 
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where x are the values of the states at the final time, and a is a set of Lagrange multipliers 
used to ensure that the states are continuous at the initial and final times. 

Next, to simplify notation, we define a Hamiltonian as 

H = L + \ T f + n T g (6) 


and define xj = x(tf) and xo = x(£o) To satisfy the first-order necessary conditions for a 
local m inimu m , we follow the development in [1] and take the first variation of J, allowing 
variations in the states, state rates, controls, Lagrange multipliers, slack variables, and 
final time, yielding 


6J* = ^^ + v T ^ + L + \ T (f-x) + ti T (g + k?)^ dt f + 6v T '& 

+ (&7 + vT w) ^ dx((o) + 6a{i - 1)! ‘- 

+ a{dx - dx)\\i 4- £' [*** + - A T Sx + SX T (f - x) 

4- 6n T (g 4 k 2 ) + 26fi T k6kj dt 


(7) 


Given that the variations are continuous, dx = dx. We then choose Sot = dX, integrate 
X T Sx by parts and expand the total differentials at the end points. 


dx{to) — Sx(to) 

dx(tf) = 6x(tf) + x(tf)dtf 

dX(t f ) = 6X(t f ) + X(t f )dt f 


Eq. (7) then becomes 


d<f> 


dv 


"'"*(£ *-';s 


) 


6x(to) 


+ dtf 


K 


dxj dx / J 


4- A T (x — x) — A 7 x 4- n 1 k 


T U 2 


n f 


+ dt/ + i/ T ^- + 4- Sv T V + SX T (x - x)\l f 0 - X T Sx\\ f 0 


*/ 

+ f 1 [ H 6x + 4- X T Sx + SX T (f - x) 

Jt 0 [dx du 

4- Sfi T (g -I- k 2 ) 4- 2 fi T kSk] dt 


( 8 ) 


(9) 
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Defining subscripts on H to denote partial derivatives and rearranging terms gives: 

6J ' = (!* + " T fr + ff ) ( , + * /iT + - A ) t/ 

+ (1^ + " T i^ ■ AT ) t/ &(t/) + (I£ + + A ^) fe(to) 


4- dtf\ T (x — x) + dt/^k^ltj 4- 8v T ^! 4- 6A t (x — x)|^ 4- f [H u 8u 

Jt 0 

4- H x 8x 4- X t 6x 4- 8X T (f — x) + 8fi T (g 4- k 2 ) 4- 2/z T A:<5A:j dt 


Now, defining A as 


;r _ 


d<j> T 
— v 


dV 


dx 0 


dx 0 

d(j> T dV 
+ u T 


t = t o 

t = t f 


{ dxf Ox/ 

and integrating the terms X T 6x and SX T x by parts, Eq. (10) becomes: 

6f = dt } + u T ^ + H^j 4- X t 6x\Z - 6X T x\ t t f 0 

■ rp -1 £y 

4- ^A — A^ 6x 4- 8v T ty + dtf |± T (A — A) J 
J to 

— [(x — x) t 6A]^ 4- dtfH T k 2 \ tf — dtf ^A r (x — x)j 

4 - J ' \h u 8u + HJx - X t 8x + 8X T f 4 - 8X T x 
+ 8n T (g 4- k 2 ) + 2 fx T k8k] dt 

We will now enforce that x = x and A = A at to and t / and define 

- d<j> T d 4' 

Hs -£ + 1 ' T> 7 

then Eq. (12) becomes 

8J' = dtf (jff + H ^ 4- dtffi T k 2 \t f + 8v T 't' 4- A T 6x|^ — x T ^A|^ 

4- J 4- H x 8x — 8x t X 4- 8X T f 4- 8X T x 

4- 8fi T (g 4- A: 2 ) 4- 2fi T k8k] dt 


( 10 ) 


( 11 ) 


( 12 ) 


(13) 


( 14 ) 
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The time interval is broken up into N not necessarily equal length time elements, Atj, 
such that the time at each element boundary U is calculated as: 

= < -° , (15) 

ti = t,— i + Ati i = 2, . . . , N -+■ 1 

and we define the states and controls at these nodes to be 

Xi=x(ii) i = l,...,N+l (16) 

Ui = u(ti) i=l,...,N+l (17) 

The time within the i th element, is expressed as U = U-\ + rAt* where 0 < r < 1 and 

r = hzJh=L (18) 

Ati 

so that dt = Atidr. Substituting these relationships into Eq. (14) gives: 

6f — dtf (fl 4 - + dt f n T k 2 \tj + 6u T V + X T 6x\\ J 0 - x T 6\(* 

+ £>, £ [*„*« + «•„«*, + S -^x, - ( 19 > 

4- 6 Af /* + SfjtJ (< gi + fc?) + 2 /if dr 

where subscript i refers either to quantities within the i th element or to functions evaluated 
using quantities within the i th element. 

Higher-Order Shape Functions 

FVom [2], we define C° shape functions for the variations of state equations’ Lagrange 
multipliers (or costates) and states in each element in terms of nodal values (superscript ') 
and internal values (superscript - ). 


SXi = 6Xi(l - r) + <5A i+ ir + (1 - T)Tj8,-(r)tfAy 

( 20 ) 

ri6~ 1 

SXi = 6£i(l - t ) + SXi+iT + ^2 (1 - T)TPj(T)6Xij 

i= 1 

Here rib is the order of the shape function polynomial being used, with rib = 1 representing 
h-version shape functions, and the summation would be ignored. The functions 0j(r) are 
polynomials of order (J — 1) as defined in [2]. That functional form is necessitated by the 
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time derivatives in Eq. (19) and required end conditions. The time derivative expression 
looks like 

ddx- ” t_1 

= Sx'i = AtiSii = -bii + 6x i+ 1 + 22 jj(r)6xij (21) 

T 3=1 

where 

1 jW = [(1 - r)r^(r)]' 

Similar expressions hold for the costates. 

No time derivatives of the states and costates themselves are needed in Eq." (19), so 
simpler shape functions are used 


Xi t = 0 

y^gj(r)A 0 - 0 < r < 1 (22) 

i=l 

A i+ l T = 1 

where the functions aj (r) are again polynomials of order j — 1 as defined in [2]. 

For the controls, control constraint Lagrange multipliers, slack variables, and their 
variations, again no time derivatives exist in Eq. (19), so the same easier shape functions 
axe used: 


r»6 




6uj = 22 a - 

3=1 

}(r)6uij 

(23) 



nb 

nb 



6 ^ = 


6ki = 22 a j( T )6kij 

(24) 



i= i 

3=1 




nb 





Ui =22 aj(r)uij 
3=1 

0 < r < 1 

(25) 

nb 



nb 


Vi =22 


0 < T < 1 

ki = 22 a 3( T )kij 0 < r < 1 

(26) 


3=1 3=1 


Xi = i 

r 

nb 

E 

Xi 

Ctj(T)Xij 

T = 0 

0 < r < 1 

Ai 


j= 1 

£i+l 

T = 1 
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Substituting these relationships into Eq.(19) gives 

6J' — + Ajv +1 &r/v+i — XfSxi — £jv+i^Ajv+i 4- x^Ai 

. N rl ( / "6 X 

+ 6v T * + dt f (H+H) t ,+Y j { AU ^ | H Ui 

“ ^ 

Tlfc — 1 

<5ii(l - t) + 6£ i+1 T + €j(T)6Xij 


+ H Xi 

+ J?\ 


i=i 

tit — 1 


tit— l 

<5Aj(l - r) + SXi+ir + ^ e_,(r)5Ay 

i=i 



nfc 

T 

nfc 

+ 

X] Oj (r )6fiij 




i= 1 


i=i 


+ 2 


+ 


( n t \ T / n 6 \ n fc | 

XJ Oj J I a i ( T )hj j J2 a i ( T ) 6 *ij ? dT ) 

J „ -6-1 lT - 
—6Xi + 6A i+ i + 7jf( r )^Ai J - 

-I 


* ri 

?/ 


N ,1 


nfc — 1 
/ 

3 - 1 

nfc — 1 


nfc 

J2<Xj ( T )zij dr 

3 = 1 


nfc — i Ufc 

4- <5xj-|_i -|- ^ ^ 7j ^ ] Q j ) Ajj dr 

i = i j'=i 


(27) 


where 

«i(*0 = (1 - r)r/3j(r) 


Due to the orthogonality of the chosen polynomials a(r) and 7 (r), this equation 
reduces to 
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6J' — <f£//2w+i^]v+i + ^N+if> x N+i — Af&ri — x]v +1 5Ajv+i +x^5Ai 
4- 6v t V + dt f [H + if) t + £ Ati jf 1 1 if ttl ^ (Tjfliy j 
+ 6x [ (1 - r)H Xi + 6xJ +l TH Xi + 6Xf (1 - r)fi + S\J +1 Tfi 

nt, — 1 

+ X] [ e i( T ) 6i Ij H xi + tj(T)6>Zjfi] 

J = 1 


+ I 53 ( T )^i I & + 53 a i ( T ) k ij 

\j = i / i=i 


+ 2 1 53 I 1 53 a i( T )*y I | dr 


+ 53 

< V i=2 


+ I ^A it i - ^^.iAi.i + y^6xi,j-iAjj 

i V J=2 
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Rearranging terms gives 


6J' =6Xl 

N 


dr 


-Xi + Xi t i - A*1 f (1 - r)/i 

Jo 

+ y^6\J -Xi-1,1 - Ati-1 [ rfi-1 dr + Xi'i - At* [ (1 -r)/idr 

^ L Jo Jo J 


+*A r 


N+l 


-£jv,i - Atjv / r/jv dr + xjv+i 

Jo 


rsxj. 

[- -Zi.i+1 + Afi f tj(T)fi dr 



=i 

L Jo 

■ 



( . cl 



-\-6x[ < 

— Ai ( i — Atj J (1 — t) 

L *i 

j=i 

dr > 

J 


N 

r 

r 1 

n 6 

+ XX < 

SO 

<1 

1 

•*** 

l-< 

1 

tO-r) 

L *< +S a i( T )^5/*€ 

»=2 

l 

Jo 

:?=i 


dr 


4- Ai-i,i — At 






i=i 


dr 


TV n 6 -l 

r ci 



+EE H «' 

Aij+i + Ati / ej(r) 


dr ► 

i=l J = 1 

J 0 

i =1 

> 


+<5x^ + i ^ -Ajv+i + A^r,i - At at / r 


r 


n* 


i=i 




dr 


+dt//t^ + ifc^ + i 4- <5 jc 7 W 

JV n 6 


\ k r 1 

+dt f (H + H) i +^J26u ij At i J H Ui CLj{T).dT 


N n (, -1 

/ a j( T ) 

< ^ • / ° 

iV nt .1 

+EE^ Ai ‘ / 2Q i< r ) 

i j=i • / ° 


«6 


i= i 


dr 


nt 

T 

nb 

X^ a i( T )£v 


y~] Q j( T )kj 

_i=i 


.j =i 


dr 


(29) 


This grouping of terms shows the equations to be solved as coefficients of the various 
variational quantities. When the coefficients are set equal to zero, the variation of J' is 
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zero, approximating the first-order necessary conditions for optimal control. In addition to 
the above, the costate boundary conditions need to be enforced, as provided for previously 



d<j> 

dx 0 


It 

Xn+i ~ &r f 



Oxq 

dxf 


= 0 
= 0 


(30) 


The equations when there axe multiple phases are similar and can also be handled 
by the code. Extra boundary conditions have to be specified, and the appropriate jump 
conditions for the costates and Hamiltonian are handled automatically. 


Implementation 

The above equations axe solved using a restricted-step Newton-Raphson method, as 
implemented in a FORTRAN code. Sparse linear systems solvers from the Harwell subrou- 
tine library [3] axe used. The user needs to specify an initial guess for all of the variables 
for the Newton-Rapheson iteration. 

The symbolic-manipulation package Macsyma developed by Symbolics [4] is used to 
generate analytic partial derivatives. Macsyma is also used to generate the a(r) and e(r) 
polynomials, which come from a recursion formula involving derivatives and integrals of 
polynomials, as developed in [2]. The user specifies the order of the polynomials, and 
Macsyma generates the necessary ones. At this time all variables have the same order 
shape function, but eventually the order for each variable will be able to be independently 
specified. The order of the shape functions can be changed between runs, but only if the 
desired number is less than that specified when the Macsyma-generated code was made. 
Otherwise the necessary polynomials will not be available. 

The integrals in Eq. 29 are approximated using Gaussian integration, with the user 
selecting the number of Gauss points, which at this time is constant for all of the inte- 
grations. That number, as well as the number of elements can be adjusted between runs 
through use of a namelist file and appropriate changes to the initial guess file. 


Results 

The code has been tested on a linear, single-state, single-control system and on a 
multiple-state, single-control system with nonlinear system dynamics. 

Linear Problem The linear problem we considered is a minimum energy problem for getting 
from one position to another. The system dynamics axe 

x = x + u (31) 
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where x is the state and u is the control. The boundary conditions are 


^1 = xo — (e - 1) 

= X f - (1 - e) = 0 

and the cost function is 

J= f \u 2 dt 

Jo 2 

The analytic solution to this problem is 

x(t) = e^ 1- ^ — e* 
A(t) = 2e (1_t) 
u(t) = — 


( 32 ) 


(33) 


(34) 


This problem was run for a variety of shape function orders, numbers of Gauss points, 
and numbers of elements. The code assumes a free time problem, introducing a nonlinear 
equation ( H{tf ) — H = 0) even into linear problems. Knowing this, in all cases all of the 
linear equations were solved in a single iteration, while that last nonlinear equation con- 
verged better than quadratically. This problem is easy enough to work by hand reasonably, 
and the Jacobian matrix generated by the code matched the one we derived explicitly for 
the second-order shape function case. 

Also, since all of the errors in all of the equations in (29) were solved to within l.Oe— 10, 
the values for the costate were always just the negative of those for the control and the 
state boundary conditions were met. Therefore in comparing results, we will look at the 
initial and final values of only the control. 

Table 1: Errors in initial and final values of control 
for various values of higher-order finite-element parameters 


u(to) 

Error 

U(tf) 

Error 

Gauss 

Points 

Shape Fn. 
Order 

Number of 
Elements 

5.18e-2 

1.41e-l 

1 

1 

1 

8.58e-4 

2.33e-3 

2 

2 

1 

6.00e-6 

1.63e-5 

3 

3 

1 

2.36e-8 

6.41e-8 

4 

4 

1 

5.92e-ll 

1.61e-10 

5 

5 

1 

1.23e-2 

3.34e-2 

2 

1 

2 

5.12e-5 

1.39e-4 

2 

2 

2 

9.12e-8 

2.48e-7 

3 

3 

2 

9.01e-ll 

2.45e-10 

4 

4 

2 

5.35e-14 

1.50e-13 

5 

5 

2 
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As one can see, the errors reduce dramatically as the order of the shape functions 
increases. Not shown is a similar decrease in error as the number of Gauss points increases, 
to a minimum when the number of Gauss points equals the order of the element. This is 
similar performance to the time-marching algorithm, as presented in our previous report. 
With further analysis we will do trade studies comparing the computational effort necessary 
to solve the problem with more elements vs. using higher-order shape functions. 


Nonlinear Problem, Next we computed some preliminary solutions to a problem with 
n onlinear system dynamics that involves maximum velocity transfer of a particle of mass 
m to a rectilinear path in a fixed time (see [5], pg. 59). The mass is acted on by a force of 
constant magnitude ma and variable heading /?(<). The states for this problem are position 
of the particle x (horizontal) and y (vertical) and the corresponding velocity components 
ti and v. The differential equations for this system are then: 


x — u 
y = v 
u — a cos (3 
v = asin/3 


(35) 


with initial conditions corresponding to a zero velocity at the origin and with terminal 
constraints that the particle is in horizontal flight at a given height (assumed here to be 
1) 

V i = u(0) = 0 

^2 = v(0) = 0 


^3 = x(0) = 0 

*4 = 2 /( 0 ) = 0 


(36) 


^5 = y(tf) -i = o 

$6 = v(tf) — o 


for an unspecified final horizontal position x(t j) and for a final horizontal velocity u(t /) 
to be maximized. The cost function is then 


J = u(tf) 


(37) 


Reference [5] gives the analytic solution in terms of the initial force heading angle, the 
final time, and the final altitude in unspecified units. These values were chosen to be 75° , 1, 
and 1 respectively. In all cases, all of the boundary conditions were met and the horizontal 
velocity was twice the horizontal position, leaving the final values of horizontal velocity, 
horizontal position costate and heading angle as the three most interesting quantities to 
look at for comparison purposes. 
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Table 2: Errors in final values of U(tf ), A „($/), and /3(tf) 
for various values of higher-order finite-element parameters 


U(tf) 

error 

Aj ,(*/) 
error 

/%) 

error 

Number of 
Elements 

Shape Fn. 
Order 

1.90e-01 

4.25e-02 

8.07e-03 

2 

1 

2.77e-02 

7.34e-02 

1.39e-02 

2 

2 

N/A 

N/A 

N/A 

3 

1 

1.89e-03 

7.79e-03 

1.49e-03 

3 

2 

6.77e-02 

7.53e-02 

1.42e-02 

4 

1 

2.80e-03 

1.55e-02 

2.95e-03 

4 

2 

5.31e-02 

1.21e-01 

2.28e-02 

5 

1 

4.77e-04 

1.52e-03 

2.90e-04 

5 

2 


Pr elimin ary results axe shown in Table 2. Missing values were for cases that did not 
converge, which is not surprising given the small number of elements being used. For this 
nonlinear problem the improvement with the higher-order element is not as dramatic. For 
4 and 5 elements all errors decrease with the higher-order elements. For 2 elements some 
errors do not decrease. This may be due to the need for a larger number of Gauss points 
when the order is increased (6 points were used). For the simplest (h-version) element, 
1 Gauss point is optimum in the sense that it is the smallest number of points which 
gives acceptable error. However, when one adds the higher-order shape functions, a larger 
number of Gauss points will always be required for nonlinear problems than for linear ones. 

The optimum number of Gauss points seems to be dependent on both the maximum 
element order and the types of nonlinearities in the problem. It is possible that there is a 
simple calculation, using some terms from the equations under consideration, which can be 
done for any given problem and which will approximately determine the optimum number 
of Gauss points for that particular problem as a function of the orders of element from 2 
up to the maximum needed. Whether this is in fact true should become clearer as more 
problems are tried and more cases of these problems are run. Recent developments by 
Hinnant (of the Computational Mechanics Branch at Langley) may also lead to savings in 
the numerical quadrature costs. 

One possible drawback of higher-order finite element schemes is the increased com- 
putational effort within each element required in implementing hp-ve rsion finite elements. 
We will ultimately determine whether this computational effort is sufficiently offset by the 
reduction in the number of time elements used and improvements in the Newton-Raphson 
convergence so as to be useful in solving optimal control problems in real time. Also, 
because certain of the element interior unknowns can be eliminated at the element level by 
solving a small set of nonlinear algebraic equations in which the nodal values are taken as 
given, the scheme may turn out to be especially powerful in a parallel computing environ- 
ment. A different processor could be assigned to each element. The number of processors, 
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strictly speaking, would not be required to be any larger than the number of sub-regions 
which are free of discontinuities. 
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